libname MODEL "W:\Curdy\SAS datasets\DATASETS FOR MODELS\revised datasets with alter IDs";
data work.CC_DD_LONG_d; set MODEL.model_with; run;
data work.CC; set MODEL.CC; run;
data work.DD; set MODEL.DD; run;
/*######################################l#####################################################################*/
/*need to combine APID and AOID into one variable*/

/*create a NUMERIC dummy for CC/DD*/
data work.CC_DD_long_d; set work.CC_DD_long;
CC=0; if Origin="District to District (USFS)" then CC=1; run;

/*create a new alterID variable  that is combination of AOID (county alter) and APID (Forest)*/
 
data work.cc_dd_long_d;set work.CC_DD_long_d;
Alter_ID = AOID;
run;

data work.cc_dd_long_d;set work.CC_DD_long_d;
   if missing(Alter_ID) then Alter_ID = coalesce(APID,.);
run;


/*FEDERAL social capital expectancy*/

/*LINEAR*/
proc mixed data=DD METHOD=ML COVTEST;
CLASS EgoID;
model T1_conf = T1_IntFreq / SOLUTION CL; 
random intercept / subject=EgoID;
run;
/*QUADRATIC*/
proc mixed data=DD METHOD=ML COVTEST;
CLASS EgoID;
model T1_conf = T1_IntFreq|T1_IntFreq / SOLUTION CL; 
random intercept / subject=EgoID;
run;


/*graphing averages over levels of T1_intfreq*/

proc mixed data=DD METHOD=ML COVTEST;
  CLASS EgoID T1_intfreq;
  model T1_conf = T1_intfreq / SOLUTION; 
  random intercept / subject=EgoID;
  lsmeans T1_intfreq / pdiff cl;
run;


/*graphing probabiliitesover levels of t1_infreq*/

/*running the GLIMMIX on T1_conf_3p - 1 cofidence variable is T1_conf_3p which means 1=confidence level 3 or greater (versus 0= missing, level 0, level 1, or level 2). creating the LoO for to get the fixed effect probability
note - this code is erroring out - need ot figure out why*/
/* Step 1: Run the PROC GLIMMIX procedure */
proc glimmix data=DD;
    class T1_IntFreq (ref='0') EGOID;
    model T1_conf_3p = T1_IntFreq / dist=bin link=logit solution;
    random EGOID / solution;
    /* Output predicted probabilities and confidence limits for each observation */
    output out=morphpred pred(ilink)=pred_prob residual=resid ucl(ilink)=up lcl(ilink)=low;
    /* Calculate least-squares means (average probability) for each level of T1_IntFreq with confidence limits */
    lsmeans T1_IntFreq / ilink cl adjust=tukey;
    /* Save the output datasets */
    ods output LSMeans=T1_conf_3pmean_probs; /* LSMeans for T1_IntFreq */
    ods output solutionR=RE_sig; /* Random effect solutions */
    ods output ParameterEstimates=LO_est; /* Parameter estimates */
    ods output diffs=diff_level; /* Differences of LSMeans */
run;

/* Step 2: Process diff_level dataset to extract specific pairs and differences */
data stats;
    retain Pair PP_diff DF Pr_diff;
    set diff_level;
    if _N_ in (2, 3, 4, 7, 8, 9, 11, 12, 14, 15) then delete;
    if _N_=1 then Pair="12";
    else if _N_=6 then Pair="23";
    else if _N_=10 then Pair="34";
    else if _N_=13 then Pair="45";
    else if _N_=5 then Pair="01";
    PP_diff = round(Estimate, .001);
    if Probt < .0001 then Pr_diff = .0001;
    else Pr_diff = round(Probt, .0001);
    keep Pair PP_diff DF Pr_diff;
run;

/* Step 3: Creating the margins1 dataset from LSMeans data */
data margins1;
    retain Marg_Level Marg_eff Pr_Marg_eff;
    set T1_conf_3pmean_probs;
    Marg_Level = T1_IntFreq;
    Marg_eff = round(Estimate, .0001);
    if Probt < .0001 then Pr_Marg_eff = .0001;
    else Pr_Marg_eff = round(Probt, .0001);
    keep Marg_Level Marg_eff Pr_Marg_eff;
run;

/* Step 4: Process the RE_sig dataset to retain random effects estimates */
data Re_sig1;
    retain RE RE_est Pr_RE;
    set RE_sig;
    RE = Effect;
    RE_est = Estimate;
    Pr_RE = Probt;
    keep RE RE_est Pr_RE;
run;

/* Step 5: Process the LO_est dataset to retain fixed effect parameter estimates */
data LO_est1;
    retain order LO_level LoO Pr_LO;
    set LO_est;
    if Effect = "Intercept" then LO_level = "Intercept";
    else LO_level = put(T1_IntFreq, 1.);
    if LO_level = "Intercept" then order = 1;
    else order = input(LO_level, best.);
    LoO = round(Estimate, .001);
    if Probt < .0001 then Pr_LO = .0001;
    else Pr_LO = round(Probt, .0001);
    keep order LO_level LoO Pr_LO;
run;

/* Step 6: Calculate random effect estimate LoO as probabilities (0-100) for each level of T1_IntFreq with confidence limits */
data T1_conf_3pLoO_prob;
    set LO_est1;
    if LO_level = "Intercept" then LO_level = "Intercept";
    /* Ensure StdErr is correctly defined or calculated; for this example, assuming it's not available */
    /* Use an arbitrary value if StdErr is not available; otherwise, correct this part based on your data */
    StdErr = 0.1; /* Adjust or remove based on available data */
    LoO_prob = round((exp(LoO) / (1 + exp(LoO))) * 100, .1); /* Convert to percentage */
    LoO_lcl_prob = round((exp(LoO - 1.96 * StdErr) / (1 + exp(LoO - 1.96 * StdErr))) * 100, .1); /* Lower confidence limit */
    LoO_ucl_prob = round((exp(LoO + 1.96 * StdErr) / (1 + exp(LoO + 1.96 * StdErr))) * 100, .1); /* Upper confidence limit */
    keep LO_level LoO_prob LoO_lcl_prob LoO_ucl_prob Pr_LO;
run;

/* Step 7: Combine all processed datasets into DD_T1_conf_3p_margins */
data DD_T1_conf_3p_margins;
    retain DV Diffs Pair DF logit_diff Pr_diff Marg_Level Marg_eff Pr_Marg_eff RE RE_est Pr_RE Level Prob Pr;
    /* Initialize with the stats dataset */
    merge stats(in=a) margins1(in=b) Re_sig1(in=c) LO_est1(in=d);
    by order;
    if a and b and c and d;
    /* Compute Marg_eff and Prob from logit values */
    if Marg_eff ^= . then Marg_eff = round(exp(Marg_eff) / (1 + exp(Marg_eff)), .001);
    if LoO ^= . then Prob = round(exp(LoO) / (1 + exp(LoO)), .001);
    Level = LO_level;
    Pr = Pr_LO;
    logit_diff = PP_diff;
    drop PP_diff LO_level LoO Pr_LO;
run;

/* Step 8: Append additional datasets to the stats dataset */
proc append base=stats data=margins1 force; run;
proc append base=stats data=Re_sig1 force; run;
proc append base=stats data=LO_est1 force; run;

/************************County Social Capital Expectancy*/
/*testing for curvilinearity*/


proc mixed data=CC METHOD=ML COVTEST;
CLASS EgoID AOID;
model T1_nfrnd = T1_IntFreq / SOLUTION cl; 
random intercept / subject=EgoID;
random intercept / subject=AOID;
run;
/*QUADRATIC*/
proc mixed data=CC METHOD=ML COVTEST;
CLASS EgoID AOID;
model T1_nfrnd = T1_IntFreq|T1_IntFreq / SOLUTION cl; 
random intercept / subject=EgoID;
random intercept / subject=AOID;
run;

/*  /*graphiing averages over levels of intfreq*/
proc mixed data=CC METHOD=ML COVTEST;
  CLASS EgoID T1_intfreq;
  model T1_nfrnd = T1_intfreq / SOLUTION; 
  random intercept / subject=EgoID;
  lsmeans T1_intfreq / pdiff cl;
run;



/**graphing probabilities**/
/*T1_nfrnd_4p   =1 if T1_nfrnd is equal to 4 or 5, else T1_nfrnd_4p=0 .



/* Step 1: Run the PROC GLIMMIX procedure */

proc glimmix data=CC;
    class T1_IntFreq (ref='0');
    model T1_nfrnd_4p = T1_IntFreq / dist=bin link=logit solution;
    random EGOID AOID/ solution;
    /* Output predicted probabilities and confidence limits for each observation */
    output out=morphpred pred(ilink)=pred_prob residual=resid ucl(ilink)=up lcl(ilink)=low;
    /* Calculate least-squares means (average probability) for each level of T1_IntFreq with confidence limits */
    lsmeans T1_IntFreq / ilink cl adjust=tukey;
    /* Save the output datasets */
    ods output LSMeans=mean_probs; /* LSMeans for T1_IntFreq */
    ods output solutionR=RE_sig; /* Random effect solutions */
    ods output ParameterEstimates=LO_est; /* Parameter estimates */
    ods output diffs=diff_level; /* Differences of LSMeans */
run;

/* Step 2: Process diff_level dataset to extract specific pairs and differences */
data stats;
    retain Pair PP_diff DF Pr_diff;
    set diff_level;
    if _N_ in (2, 3, 4, 7, 8, 9, 11, 12, 14, 15) then delete;
    if _N_=1 then Pair="12";
    else if _N_=6 then Pair="23";
    else if _N_=10 then Pair="34";
    else if _N_=13 then Pair="45";
    else if _N_=5 then Pair="01";
    PP_diff = round(Estimate, .001);
    if Probt < .0001 then Pr_diff = .0001;
    else Pr_diff = round(Probt, .0001);
    keep Pair PP_diff DF Pr_diff;
run;

/* Step 3: Creating the margins1 dataset from LSMeans data */
data margins1;
    retain Marg_Level Marg_eff Pr_Marg_eff;
    set t1_nfrnd_4pmean_probs;
    Marg_Level = T1_IntFreq;
    Marg_eff = round(Estimate, .0001);
    if Probt < .0001 then Pr_Marg_eff = .0001;
    else Pr_Marg_eff = round(Probt, .0001);
    keep Marg_Level Marg_eff Pr_Marg_eff;
run;

/* Step 4: Process the RE_sig dataset to retain random effects estimates */
data Re_sig1;
    retain RE RE_est Pr_RE;
    set RE_sig;
    RE = Effect;
    RE_est = Estimate;
    Pr_RE = Probt;
    keep RE RE_est Pr_RE;
run;

/* Step 5: Process the LO_est dataset to retain fixed effect parameter estimates */
data LO_est1;
    retain order LO_level LoO Pr_LO;
    set LO_est;
    if Effect = "Intercept" then LO_level = "Intercept";
    else LO_level = put(T1_IntFreq, 1.);
    if LO_level = "Intercept" then order = 1;
    else order = input(LO_level, best.);
    LoO = round(Estimate, .001);
    if Probt < .0001 then Pr_LO = .0001;
    else Pr_LO = round(Probt, .0001);
    keep order LO_level LoO Pr_LO;
run;

/* Step 6: Calculate random effect estimate LoO as probabilities (0-100) for each level of T1_IntFreq with confidence limits */
data t1_nfrnd_4pLoO_prob;
    set LO_est1;
    if LO_level = "Intercept" then LO_level = "Intercept";
    /* Ensure StdErr is correctly defined or calculated; for this example, assuming it's not available */
    /* Use an arbitrary value if StdErr is not available; otherwise, correct this part based on your data */
    StdErr = 0.1; /* Adjust or remove based on available data */
    LoO_prob = round((exp(LoO) / (1 + exp(LoO))) * 100, .1); /* Convert to percentage */
    LoO_lcl_prob = round((exp(LoO - 1.96 * StdErr) / (1 + exp(LoO - 1.96 * StdErr))) * 100, .1); /* Lower confidence limit */
    LoO_ucl_prob = round((exp(LoO + 1.96 * StdErr) / (1 + exp(LoO + 1.96 * StdErr))) * 100, .1); /* Upper confidence limit */
    keep LO_level LoO_prob LoO_lcl_prob LoO_ucl_prob Pr_LO;
run;

/* Step 7: Combine all processed datasets into CC_t1_nfrnd_4p_margins */
data CC_t1_nfrnd_4p_margins;
    retain DV Diffs Pair DF logit_diff Pr_diff Marg_Level Marg_eff Pr_Marg_eff RE RE_est Pr_RE Level Prob Pr;
    /* Initialize with the stats dataset */
    merge stats(in=a) margins1(in=b) Re_sig1(in=c) LO_est1(in=d);
    by order;
    if a and b and c and d;
    /* Compute Marg_eff and Prob from logit values */
    if Marg_eff ^= . then Marg_eff = round(exp(Marg_eff) / (1 + exp(Marg_eff)), .001);
    if LoO ^= . then Prob = round(exp(LoO) / (1 + exp(LoO)), .001);
    Level = LO_level;
    Pr = Pr_LO;
    logit_diff = PP_diff;
    drop PP_diff LO_level LoO Pr_LO;
run;

/* Step 8: Append additional datasets to the stats dataset */
proc append base=stats data=margins1 force; run;
proc append base=stats data=Re_sig1 force; run;
proc append base=stats data=LO_est1 force; run;


/***Disaster Interaction Frequency*/


proc mixed data=CC_DD_long_d METHOD=ML COVTEST;
CLASS EgoID CC;
model T2_IntFreq = T1_IntFreq CC/ SOLUTION CL; 
random intercept / subject=EgoID;
run;
/*QUADRATIC*/
proc mixed data=CC_DD_long_d METHOD=ML COVTEST;
CLASS EgoID CC;
model T2_IntFreq = T1_IntFreq|T1_IntFreq CC/ SOLUTION CL; 
random intercept / subject=EgoID;
run;

/**looking at both alter and ego DI***/

proc mixed data=CC_DD_long_d  METHOD=ML COVTEST;
CLASS EgoID Alter_ID CC;
model T2_IntFreq = T1_IntFreq|T1_IntFreq CC / SOLUTION cl; 
random intercept / subject=EgoID;
random intercept / subject=Alter_ID;
run;

proc mixed data=CC_DD_long_d  METHOD=ML COVTEST;
CLASS AlterID CC;
model T2_IntFreq = T1_IntFreq|T1_IntFreq CC / SOLUTION cl; 
random intercept / subject=Alter_ID;
run;

/***graphing means over levels of t1int_freq**/

proc mixed data=CC_DD_long_d METHOD=ML COVTEST;
  CLASS EgoID T1_intfreq;
  model T2_intfreq = T1_intfreq / SOLUTION; 
  random intercept / subject=EgoID;
  lsmeans T1_intfreq / pdiff cl;
run;



/*GLIMIX*/
/*graphiing averages over levels of intfreq*/
/*NOTE- GLIMIX model would not converge alter_ID and egoID as random effects - confirmed the results are consistent if control for Alter_ID as Ego_ID*/

proc glimmix data=CC_DD_long_d;
    class T1_IntFreq (ref='0');
    model t2_intfreq_3p = T1_IntFreq / dist=bin link=logit solution;
    random EGOID / solution;
    /* Output predicted probabilities and confidence limits for each observation */
    output out=morphpred pred(ilink)=pred_prob residual=resid ucl(ilink)=up lcl(ilink)=low;
    /* Calculate least-squares means (average probability) for each level of T1_IntFreq with confidence limits */
    lsmeans T1_IntFreq / ilink cl adjust=tukey;
    /* Save the output datasets */
    ods output LSMeans=t2_intfreq_3pmean_probs; /* LSMeans for T1_IntFreq */
    ods output solutionR=RE_sig; /* Random effect solutions */
    ods output ParameterEstimates=LO_est; /* Parameter estimates */
    ods output diffs=diff_level; /* Differences of LSMeans */
run;

/* Step 2: Process diff_level dataset to extract specific pairs and differences */
data stats;
    retain Pair PP_diff DF Pr_diff;
    set diff_level;
    if _N_ in (2, 3, 4, 7, 8, 9, 11, 12, 14, 15) then delete;
    if _N_=1 then Pair="12";
    else if _N_=6 then Pair="23";
    else if _N_=10 then Pair="34";
    else if _N_=13 then Pair="45";
    else if _N_=5 then Pair="01";
    PP_diff = round(Estimate, .001);
    if Probt < .0001 then Pr_diff = .0001;
    else Pr_diff = round(Probt, .0001);
    keep Pair PP_diff DF Pr_diff;
run;

/* Step 3: Creating the margins1 dataset from LSMeans data */
data margins1;
    retain Marg_Level Marg_eff Pr_Marg_eff;
    set t2_intfreq_3pmean_probs;
    Marg_Level = T1_IntFreq;
    Marg_eff = round(Estimate, .0001);
    if Probt < .0001 then Pr_Marg_eff = .0001;
    else Pr_Marg_eff = round(Probt, .0001);
    keep Marg_Level Marg_eff Pr_Marg_eff;
run;

/* Step 4: Process the RE_sig dataset to retain random effects estimates */
data Re_sig1;
    retain RE RE_est Pr_RE;
    set RE_sig;
    RE = Effect;
    RE_est = Estimate;
    Pr_RE = Probt;
    keep RE RE_est Pr_RE;
run;

/* Step 5: Process the LO_est dataset to retain fixed effect parameter estimates */
data LO_est1;
    retain order LO_level LoO Pr_LO;
    set LO_est;
    if Effect = "Intercept" then LO_level = "Intercept";
    else LO_level = put(T1_IntFreq, 1.);
    if LO_level = "Intercept" then order = 1;
    else order = input(LO_level, best.);
    LoO = round(Estimate, .001);
    if Probt < .0001 then Pr_LO = .0001;
    else Pr_LO = round(Probt, .0001);
    keep order LO_level LoO Pr_LO;
run;

/* Step 6: Calculate random effect estimate LoO as probabilities (0-100) for each level of T1_IntFreq with confidence limits */
data t2_intfreq_3pLoO_prob;
    set LO_est1;
    if LO_level = "Intercept" then LO_level = "Intercept";
    /* Ensure StdErr is correctly defined or calculated; for this example, assuming it's not available */
    /* Use an arbitrary value if StdErr is not available; otherwise, correct this part based on your data */
    StdErr = 0.1; /* Adjust or remove based on available data */
    LoO_prob = round((exp(LoO) / (1 + exp(LoO))) * 100, .1); /* Convert to percentage */
    LoO_lcl_prob = round((exp(LoO - 1.96 * StdErr) / (1 + exp(LoO - 1.96 * StdErr))) * 100, .1); /* Lower confidence limit */
    LoO_ucl_prob = round((exp(LoO + 1.96 * StdErr) / (1 + exp(LoO + 1.96 * StdErr))) * 100, .1); /* Upper confidence limit */
    keep LO_level LoO_prob LoO_lcl_prob LoO_ucl_prob Pr_LO;
run;

/* Step 7: Combine all processed datasets into DD_t2_intfreq_3p_margins */
data DD_t2_intfreq_3p_margins;
    retain DV Diffs Pair DF logit_diff Pr_diff Marg_Level Marg_eff Pr_Marg_eff RE RE_est Pr_RE Level Prob Pr;
    /* Initialize with the stats dataset */
    merge stats(in=a) margins1(in=b) Re_sig1(in=c) LO_est1(in=d);
    by order;
    if a and b and c and d;
    /* Compute Marg_eff and Prob from logit values */
    if Marg_eff ^= . then Marg_eff = round(exp(Marg_eff) / (1 + exp(Marg_eff)), .001);
    if LoO ^= . then Prob = round(exp(LoO) / (1 + exp(LoO)), .001);
    Level = LO_level;
    Pr = Pr_LO;
    logit_diff = PP_diff;
    drop PP_diff LO_level LoO Pr_LO;
run;

/* Step 8: Append additional datasets to the stats dataset */
proc append base=stats data=margins1 force; run;
proc append base=stats data=Re_sig1 force; run;
proc append base=stats data=LO_est1 force; run;


/****Tie Performance****/
/****testing for curvilearitye****/

proc mixed data=CC_DD_long_d METHOD=ML COVTEST;
CLASS EgoID CC;
model T2_RFI = T1_IntFreq CC/ SOLUTION; 
random intercept / subject=EgoID;
run;
/*QUADRATIC*/
proc mixed data=CC_DD_long_d METHOD=ML COVTEST;
CLASS EgoID CC;
model T2_RFI = T1_IntFreq|T1_IntFreq CC/ SOLUTION; 
random intercept / subject=EgoID;
run;

/*graphiing averages over levels of intfreq*/

proc mixed data=CC_DD_long_d METHOD=ML COVTEST;
  CLASS EgoID T1_intfreq;
  model T2_RFI = T1_intfreq / SOLUTION; 
  random intercept / subject=EgoID;
  lsmeans T1_intfreq / pdiff cl;
run;

/*graphing probabilities*/

proc glimmix data=CC_DD_long_d;
    class T1_IntFreq (ref='0');
    model t2_RFI_1p = T1_IntFreq / dist=bin link=logit solution;
    random EGOID / solution;
    /* Output predicted probabilities and confidence limits for each observation */
    output out=morphpred pred(ilink)=pred_prob residual=resid ucl(ilink)=up lcl(ilink)=low;
    /* Calculate least-squares means (average probability) for each level of T1_IntFreq with confidence limits */
    lsmeans T1_IntFreq / ilink cl adjust=tukey;
    /* Save the output datasets */
    ods output LSMeans=t2_RFI_1pmean_probs; /* LSMeans for T1_IntFreq */
    ods output solutionR=RE_sig; /* Random effect solutions */
    ods output ParameterEstimates=LO_est; /* Parameter estimates */
    ods output diffs=diff_level; /* Differences of LSMeans */
run;

/* Step 2: Process diff_level dataset to extract specific pairs and differences */
data stats;
    retain Pair PP_diff DF Pr_diff;
    set diff_level;
    if _N_ in (2, 3, 4, 7, 8, 9, 11, 12, 14, 15) then delete;
    if _N_=1 then Pair="12";
    else if _N_=6 then Pair="23";
    else if _N_=10 then Pair="34";
    else if _N_=13 then Pair="45";
    else if _N_=5 then Pair="01";
    PP_diff = round(Estimate, .001);
    if Probt < .0001 then Pr_diff = .0001;
    else Pr_diff = round(Probt, .0001);
    keep Pair PP_diff DF Pr_diff;
run;

/* Step 3: Creating the margins1 dataset from LSMeans data */
data margins1;
    retain Marg_Level Marg_eff Pr_Marg_eff;
    set t2_RFI_1pmean_probs;
    Marg_Level = T1_IntFreq;
    Marg_eff = round(Estimate, .0001);
    if Probt < .0001 then Pr_Marg_eff = .0001;
    else Pr_Marg_eff = round(Probt, .0001);
    keep Marg_Level Marg_eff Pr_Marg_eff;
run;

/* Step 4: Process the RE_sig dataset to retain random effects estimates */
data Re_sig1;
    retain RE RE_est Pr_RE;
    set RE_sig;
    RE = Effect;
    RE_est = Estimate;
    Pr_RE = Probt;
    keep RE RE_est Pr_RE;
run;

/* Step 5: Process the LO_est dataset to retain fixed effect parameter estimates */
data LO_est1;
    retain order LO_level LoO Pr_LO;
    set LO_est;
    if Effect = "Intercept" then LO_level = "Intercept";
    else LO_level = put(T1_IntFreq, 1.);
    if LO_level = "Intercept" then order = 1;
    else order = input(LO_level, best.);
    LoO = round(Estimate, .001);
    if Probt < .0001 then Pr_LO = .0001;
    else Pr_LO = round(Probt, .0001);
    keep order LO_level LoO Pr_LO;
run;

/* Step 6: Calculate random effect estimate LoO as probabilities (0-100) for each level of T1_IntFreq with confidence limits */
data t2_RFI_1pLoO_prob;
    set LO_est1;
    if LO_level = "Intercept" then LO_level = "Intercept";
    /* Ensure StdErr is correctly defined or calculated; for this example, assuming it's not available */
    /* Use an arbitrary value if StdErr is not available; otherwise, correct this part based on your data */
    StdErr = 0.1; /* Adjust or remove based on available data */
    LoO_prob = round((exp(LoO) / (1 + exp(LoO))) * 100, .1); /* Convert to percentage */
    LoO_lcl_prob = round((exp(LoO - 1.96 * StdErr) / (1 + exp(LoO - 1.96 * StdErr))) * 100, .1); /* Lower confidence limit */
    LoO_ucl_prob = round((exp(LoO + 1.96 * StdErr) / (1 + exp(LoO + 1.96 * StdErr))) * 100, .1); /* Upper confidence limit */
    keep LO_level LoO_prob LoO_lcl_prob LoO_ucl_prob Pr_LO;
run;

/* Step 7: Combine all processed datasets into DD_t2_RFI_1p_margins */
data DD_t2_RFI_1p_margins;
    retain DV Diffs Pair DF logit_diff Pr_diff Marg_Level Marg_eff Pr_Marg_eff RE RE_est Pr_RE Level Prob Pr;
    /* Initialize with the stats dataset */
    merge stats(in=a) margins1(in=b) Re_sig1(in=c) LO_est1(in=d);
    by order;
    if a and b and c and d;
    /* Compute Marg_eff and Prob from logit values */
    if Marg_eff ^= . then Marg_eff = round(exp(Marg_eff) / (1 + exp(Marg_eff)), .001);
    if LoO ^= . then Prob = round(exp(LoO) / (1 + exp(LoO)), .001);
    Level = LO_level;
    Pr = Pr_LO;
    logit_diff = PP_diff;
    drop PP_diff LO_level LoO Pr_LO;
run;

/* Step 8: Append additional datasets to the stats dataset */
proc append base=stats data=margins1 force; run;
proc append base=stats data=Re_sig1 force; run;
proc append base=stats data=LO_est1 force; run;
